SUBROUTINE G_RADENE
!-----------------------------------------------------------------------
!
! Definition von G:   Strahlungsenergiegleichung (0.tes Moment der Strahlungstransportgleichung)
!

      use primvar,  only : G, X, XA, MH, MJ, MR, MUr, Gnorm
      use physco,   only : z1, z2, z3, z12, StBolz, pi, clight
      use config,   only : np
      use global,   only : tst
      use matvar,   only : Tgas, OPAplk
      use eddvar,   only : f_edd
      use geomvar,  only : S_flux, S_vol, S_volA
      use zvar,     only : rho0
      use advecvar, only : J_adv


      implicit none

      integer       :: i
      logical, save :: init = .true.

      if (init) then
         write(66,"(a)") "G_radene.f90:   rechne voll implizit;   Innenrand: J=S ;  Aussenrand: J=S"
         init = .false.
      end if


!-----------------------------------------------------------------------
!    Randbedingungen
!-----------------------------------------------------------------------

   ! innere Pseudozellen: i=1 & i=2
      G(MJ,1)        =  X(MJ,1) - StBolz/pi * Tgas(1)**4
      Gnorm(MJ,1)    =  max( X(MJ,1), StBolz/pi * Tgas(1)**4 )
      G(MJ,2)        =  X(MJ,2) - StBolz/pi * Tgas(2)**4
      Gnorm(MJ,2)    =  max( X(MJ,2), StBolz/pi * Tgas(2)**4 )

   ! aeussere Pseudozellen + ueberzaehlige skalare Zelle: i=np, i=np-1, i=np-2
      G(MJ,np)       =  X(MJ,np) - StBolz/pi * Tgas(np)**4
      Gnorm(MJ,np)   =  max( X(MJ,np), StBolz/pi * Tgas(np)**4 )
      G(MJ,np-1)     =  X(MJ,np-1) - StBolz/pi * Tgas(np-1)**4
      Gnorm(MJ,np-1) =  max( X(MJ,np-1), StBolz/pi * Tgas(np-1)**4 )
      G(MJ,np-2)     =  X(MJ,np-2) - StBolz/pi * Tgas(np-2)**4
      Gnorm(MJ,np-2) =  max( X(MJ,np-2), StBolz/pi * Tgas(np-2)**4 )


!-----------------------------------------------------------------------
!    Restlicher Bereich
!-----------------------------------------------------------------------

      do i=3,np-3

      ! Diese version hier ist jetzt voll implizit angeschrieben
      ! - jedoch mit ausnahme des advektionsterms wo das keinen rechten sinn macht
      !  Strahlungsdruck  "K" = f_edd(i) * X(MJ,i)
      !  source function  "S" =  StBolz/pi * Tgas(i)**4
      !         by the way   StBolz/pi = 1.80469e-05

         G(MJ,i) = S_vol(i) * X(MJ,i) - S_volA(i) * XA(MJ,i)                                                         & ! temporal diff.
                 + S_flux(i+1) * J_adv(i+1) - S_flux(i) * J_adv(i)                                                   & ! advection
                 + clight * z2*pi * ( X(MR,i+1) * X(MH,i+1) - X(MR,i) * X(MH,i) ) * tst                              & ! div(H) term - radial
                 + f_edd(i) * X(MJ,i) * z2*pi * ( X(MR,i+1) * X(MUr,i+1) - X(MR,i) * X(MUr,i) ) * tst                & ! K div(u) term
                 - z12 * X(MJ,i) * ( z3*f_edd(i) - z1 ) * S_vol(i) * tst *                                           & ! (3K-J)
                      ( X(MR,i) * X(MUr,i) + X(MR,i+1) * X(MUr,i+1) ) / ( X(MR,i)**2 + X(MR,i+1)**2 )                & !        * u/r
                 + clight * OPAplk(i) * rho0(i) * ( X(MJ,i) - StBolz/pi * Tgas(i)**4 ) * S_vol(i) * tst                ! (J-S) Term

         Gnorm(MJ,i) = max( &
              abs( S_vol(i) * X(MJ,i) ),                                                                             & ! cell content
              abs( S_flux(i+1) * J_adv(i+1) ), abs( S_flux(i) * J_adv(i) ),                                          & ! advection
              abs( clight * z2*pi * ( X(MR,i+1) * X(MH,i+1) - X(MR,i) * X(MH,i) ) * tst ),                           & ! div(H) term - radial
              abs( f_edd(i) * X(MJ,i) * z2*pi * ( X(MR,i+1) * X(MUr,i+1) - X(MR,i) * X(MUr,i) ) * tst ),             & ! K div(u) term
              abs( z12 * X(MJ,i) * ( z3*f_edd(i) - z1 ) * S_vol(i) * tst *                                           & ! (3K-J)
                      ( X(MR,i) * X(MUr,i) + X(MR,i+1) * X(MUr,i+1) ) / ( X(MR,i)**2 + X(MR,i+1)**2 ) ),             & !        * u/r
              abs( clight * OPAplk(i) * rho0(i) * max( X(MJ,i), StBolz/pi * Tgas(i)**4 ) * S_vol(i) * tst )          & ! (J-S) Term
                          )

      ! In der Strahlungsenergiegleichung stellt J-S die wesentlichste Bedingung zur Bestimmung von J dar,
      ! im groessten Teil des Sterns wird J direkt durch diese Bedingung definiert - alle dynamischen terme sind dort vernachlaessigbar.
      ! Deshalb ist es notwendig den J-S Term fuer die Normierung der Strahlungsenergiegleichung aufzuspalten.

      end do


END SUBROUTINE G_RADENE


